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Gravitational waveforms from the inspiral and ring-down stages of the binary black hole coales- 
cences can be modelled accurately by approximation/perturbation techniques in general relativity. 
Recent progress in numerical relativity has enabled us to model also the non-perturbative merger 
phase of the binary black-hole coalescence problem. This enables us to coherently search for all three 
stages of the coalescence of non-spinning binary black holes using a single template bank. Taking 
our motivation from these results, we propose a family of template waveforms which can model the 
inspiral, merger, and ring-down stages of the coalescence of non-spinning binary black holes that 
follow quasi-circular inspiral. This two-dimensional template family is explicitly parametrized by 
the physical parameters of the binary. We show that the template family is not only effectual in 
detecting the signals from black hole coalescences, but also faithful in estimating the parameters of 
the binary. We compare the sensitivity of a search (in the context of different ground-based inter- 
ferometers) using all three stages of the black hole coalescence with other template-based searches 
which look for individual stages separately. We find that the proposed search is significantly more 
sensitive than other template-based searches for a substantial mass-range, potentially bringing about 
remarkable improvement in the event-rate of ground-based interferometers. As part of this work, 
we also prescribe a general procedure to construct interpolated template banks using non-spinning 
black hole waveforms produced by numerical relativity. 

PACS numbers: 



I. INTRODUCTION 

A network of ground based gravitational-wave (GW) 
detectors (LIGO, Virgo, GEO 600, TAMA) is currently 
collecting data, which a world-wide scientific collabora- 
tion is involved in analyzing. Among the most promis- 
ing sources detectable by these observatories are coa- 
lescing compact binaries consisting of black holes (BHs) 
and/or neutron stars spiraling toward each other as 
they lose orbital energy and angular momentum through 
gravitational- wave emission. The gravitational- wave sig- 
nal from coalescing binaries is conventionally split into 
three parts: inspiral, merger and ring down. In the 
first stage, the two compact objects, usually treated as 
point masses, move in quasi-circular orbits (eccentric- 
ity, if present initially, is quickly radiated away). This 
part of the waveform is described very well by the post- 
Newtonian (PN) approximation of general relativity. In 
this approximation the Einstein equations are solved in 
the near zone (which contains the source) using an ex- 
pansion in terms of the (small) velocity of the point 



masses. In the far zone, the vacuum equations are solved 
assuming weak gravitational fields, and these two solu- 
tions are matched in the intermediate region 0, H, Q . 

The PN approximation breaks down as the two com- 
pact objects approach the ultra-relativistic regime and 
eventually merge with each other. Although various re- 
summation methods, such as Pade [i| and effective-one- 
body (EOB) approaches [|[, have been developed to ex- 
tend the validity of the PN approximation, unambiguous 
waveforms in the merger stage must be calculated nu- 
merically in full general relativity. Recent breakthroughs 
in numerical relativity have allowed many groups 

@, 0, 1, S EI 03, E E| to evolve BH binaries fully nu- 
merically for the last several orbits through the plunge to 
single BH formation. The field is now rapidly developing 
the capability to routinely evolve generic black-hole bi- 
nary configurations in the comparable-mass regime, and 
to accurately extract the gravitational-wave signal. Im- 
portant milestones include simulations of unequal-mass 
binaries and calculations of the gravitational recoil ef- 
fect and the evolution of black-hole binaries with spin 
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0, HE EE El IS El HE HI, HE EE Hi . 

Comparisons with post-Newtonian results are essen- 
tial for data analysis efforts, and several groups have 
published results showing good agreement of various as- 
pects of non-spinning simulations with post-Newtonian 
predictions (see e.g. [2^, [2^, H3, [H, H^, [3^, HH, H^l ) , and 
first results for certain configurations with spin have also 
become available [IE S|- In order to overcome phase 
inaccuracies in long evolutions, significant progress has 
been made by the Caltech-Corncll group using spectral 
codes HE Hil i an d by the Jena group using higher (sixth) 
order finite differencing [37j ■ Methods to reduce the ec- 
centricity to around 10 -3 (so far only for equal-mass 
binaries) have been presented by the Caltech-Cornell 
group [36], and the Jena group [38[ (using initial pa- 
rameters from PN solutions that take into account ra- 
diation reaction). Current numerical waveforms can be 
generated for the last (< 10) orbits, and these wave- 
forms can be joined continuously with analytic PN fnspi- 
ral waveforms to obtain one full signal. This was done 
in HE HE HE HI- Indeed, there are no fundamental 
obstructions to generating the whole waveform, includ- 
ing long inspiral over hundreds of orbits, by solving the 
full Einstein equations numerically. But, not only would 
this be computationally prohibitive with current meth- 
ods, it is also unnecessary: the PN formalism is known 
to work very well in the weak-field regime (when the BHs 
are well-separated), and is a low-cost and perfectly ade- 
quate substitute to fully general relativistic solutions in 
that regime. 

The numerically generated part of the gravitational- 
wave signal from coalescing binaries also includes the 
final stage of the coalescence, when a single perturbed 
black hole is formed and it rapidly loses its deviations 
from a Kerr black hole via gravitational waves. This part 
of the signal can be decomposed as a superposition of 
exponentially damped modes, and is called quasi-normal 
mode 'ring down', by analogy with the vibrations of a 
bell. The detectable part of the ring down is rather short 
and only a few modes (if not only the dominant one) are 
expected to be important/detectable by initial ground- 
based observatories. This will not be true, however, for 
the advanced detectors and certainly it is not the 
case for LISA, the planned space-borne gravitational- 
wave observatory. Indeed, the majority of the signal- 
to-noise ratio (SNR) comes from the quasi-normal mode 
ringing of binary systems with a total mass above a few 
10 6 M Q [4l|. For LISA, and also perhaps for the next 
generation of ground based detectors, it will be possible 
to detect several quasi-normal modes and test the 'no 
hair' theorem, according to which all modes are functions 
of a BH's mass and spin [U HE El • 

Joining analytically modeled inspiral with numerically 
generated merger and ring down allows us to produce the 
complete gravitational- wave signal from coalescing bina- 
ries, and to use it in the analysis of detector data. There 
are several benefits to using the whole signal in searches. 
The most obvious one is the increase in SNR in a fully co- 



herent matched filtering search [IE EE EE 13 ■ Increase 
in SNR implies increase in the event rate and improve- 
ment in the parameter estimation. Including the inspi- 
ral. merger and ring down parts in a template waveform 
also means that the waveform has a more complex struc- 
ture. This extra complexity will also brin g ab out some 
improvement in the parameter estimation [471 ] and pos- 
sibly also a reduction in the false alarm rate in analysis 
of the data from the ground-based network of detectors. 
This is because it is in general harder for the noise to 
mimic a complex signal 1 . For LISA, the detection of 
inspiralling super-massive black holes is not a problem; 
the SNR is expected to be so large that we expect some 
signals to be visible by eye in LISA data. However, using 
the full signal for LISA data analysis is equally impor- 
tant because the full signal is essential in estimating pa- 
rameters of the binary with the required accuracy. This 
is important not only from the astrophysical point of 
view, but also because we need to subtract loud signals 
from the data in order to detect /analyze other signals. 
Imperfect signal removal due to errors in the parame- 
ter estimation will result in large residuals and will ad- 
versely affect subsequent analyses. Improved parameter 
estimation will also enable GW observations (in conjunc- 
tion with electromagnetic observations) to constrain im- 
portant cosmological parameters, most importantly the 
equation of state of dark energy [H El, EE HE HI HI . 

The numerical waveforms described above are still 
computationally expensive and cannot be used directly 
to densely cover the parameter space of the binary BHs 
that will be searched over by matched filtering tech- 
niques. A promising alternative is to use the post- 
Newtonian and numerical-relativity waveforms to con- 
struct an analytic model that sufficiently accurately 
mimics a true signal [IE H0|. In [H we have sug- 
gested a phcnomcnological family of waveforms which 
can match physical signals from non-spinning binaries 
in quasi-circular orbits with fitting factors above 99%. 
In this paper we extend this formulation to propose a 
two-parameter family of template waveforms which are 
explicitly parametrized by the physical parameters of 
the binary. We show that this two-dimensional tem- 
plate family is not only 'effectual' in detecting the sig- 
nals from binary BH coalescences, but also 'faithful' in 
estimating the parameters of the binary. This family 
of template waveforms can be used to densely cover the 
parameter space of the binary, thus avoiding the com- 
putational burden of generating numerical waveforms in 
each grid point in the parameter space. We compute the 
effectualness and faithfulness (see Section Mil for defi- 
nitions) of the template family in the context of three 
different ground-based detectors: namely, Initial LIGO, 
Virgo and Advanced LIGO. We also compare the sen- 



1 At least we expect this to happen for those binaries for which 
both the inspiral and the merger contribute significantly to SNR. 
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sitivity of a search which coherently includes all three 
(inspiral, merger and ring down) stages of the BH coa- 
lescence with other template-based searches which look 
for each stage separately. 

Our 'target signals' are constructed by matching the 
numerical-relativity waveforms to a particular family 
(TaylorTl approximant (46[) of post-Newtonian wave- 
forms, but this choice is by no means necessary. Indeed, 
we expect that more robust ways of constructing post- 
Newtonian approximants, such as the effective one-body 
approach [f| or Pade resummation approach [J] , will give 
better agreement with numerical-relativity (NR) wave- 
forms. But the purpose of the current paper is to explic- 
itly prescribe a general procedure to produce hybrid and 
phcnomcnological waveforms, and to construct interpo- 
lated template banks using parametrized waveforms. We 
show that, given the number of numerical wave cycles we 
employ, even a simple PN choice like TaylorTl leads to 
very faithful and effectual templates, and significantly in- 
creases the possible range of gravitational- wave searches. 
The use of improved PN approximants will require a 
smaller number of NR cycles, thereby further reducing 
computational cost for template construction. There are 
also other approaches for comparing analytic and numer- 
ical waveforms and for constructing hybrid waveforms 
(see, for example f2{J); it would be interesting to com- 
pare the results presented in this work with other ap- 
proaches presented in the literature. 

The paper is structured as follows. In Section lHl we 
summarize the methods of current numerical-relativity 
simulations, including a setup of the initial data that al- 
lows an unambiguous comparison with post-Newtonian 
results, and the wave extraction techniques. In Sec- 
tion IIIII we briefly outline the waveform generation us- 
ing the restricted post-Newtonian approximation. There 
wc briefly introduce the main data-analysis techniques 
and define notations that are used in the subsequent 
sections. In Section IIVI we construct a phcnomcnolog- 
ical template family parametrized only by the masses 
of the two individual black holes. First we combine re- 
stricted 3.5PN waveforms [H| with results from NR sim- 
ulations to construct 'hybrid' waveforms for the quasi- 
circular inspiral of non-spinning binaries with possibly 
unequal masses. Then, we introduce a phenomenologi- 
cal family of templates constructed in the frequency do- 
main. Initially the template family is parametrized by 
10 phcnomcnological parameters. We then find a unique 
mapping of these 10 parameters to the two physical pa- 
rameters: namely, the total mass M and the symmetric 
mass ratio r\ = M1M2/M 2 , so that the template family 
is just two-dimensional. The resulting templates have 
remarkably high fitting factors with target waveforms. 
Here we also compute the faithfulness of the templates 
and the bias in the estimation of the parameter of the 
binary. A comparison of the sensitivity of the search 
using the proposed template family with other existing 
template-based searches is also presented. Finally, we 
summarize our main results in SectionlVl Some details of 



the calculations involved are described in Appendices lAl 
and [Bj We adopt geometrical units throughout this pa- 
per: G = c = 1. 



II. NUMERICAL SIMULATIONS AND WAVE 
EXTRACTION 

Numerical simulations were performed with the BAM 
[Hi and CCATIE [H codes. Both codes evolve black- 
hole binaries using the 'moving-puncture' approach 0, 
Q. The method involves setting up initial data con- 
taining two black holes via a Brill-Linquist-like worm- 
hole construction where the additional asymptot- 
ically flat end of each wormhole is compactificd to a 
point, or 'puncture'. A coordinate singularity exists at 
the puncture, but can be stably evolved using standard 
finite-difference techniques, and is protected by causal- 
ity from adversely affecting the physically relevant ex- 
ternal spacetime. This prescription allows black holes to 
be constructed on a 3D Cartesian numerical grid with- 
out recourse to excision techniques, and also provides a 
simple way togenerate any number of moving, spinning 
black holes [HI, Hf|. Given an initial configuration of 
two black holes, the data are evolved using a conformal 
and traceless '3+1' decomposition of Einstein's equa- 
tions [13, Us L 5£ L In addition the gauge is evolved using 
the '1+log' 6j6l| and T-driver' equations [6l|,[62| and 



the coordinate singularity in the conformal factor is dealt 
with by evolving either the regular variable \ = 4'~ A 
(in BAM) or (j> = ln-0 (in CCATIE), which diverges 
'slowly' enough so as not lead to numerical instabilities. 
The standard moving puncture approach consists of all 
these techniques, and causes the 'punctures' to quickly 
assume a cylindrical asymptotics [631 ]. and allows them 
to move across the numerical grid. This method has 
been found to allow accurate, stable simulations of black 
holes over many (> 10) orbits through merger and ring 
down. 

In the initial data construction we must specify the 
masses, locations and momenta of the two black holes 
(we do not consider spinning black holes in this work). 
The mass of each black hole, Mi, is specified in terms of 
the Arnowitt-Deser-Misner (ADM) mass at each punc- 
ture. This corresponds to the mass at the other asymp- 
totically flat end which is, to a very good approximation, 
equal to irreducible mass of the apparent-horizon mass 

HEE 111 



Mi 



16tt 



(2.1) 



where Ai is the area of the apparent horizon. We assume 
that this mass is the same as the mass used in post- 
Newtonian formulas. This assumption is really expected 
to be true only in the limit where the black holes are 
infinitely far apart and stationary. As such we consider 
any error in this assumption as part of the error due 
to starting the simulation at a finite separation. The 
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important point is that a binary with horizon masses 
Mi and M 2 should be compared with a post-Newtonian 
system with the same mass parameters. This allows us 
to provide the same overall scale M = Mi + M2 for both 
numerical and post-Newtonian waveforms, and is crucial 
for comparison and matching. 

The initial momenta of the black holes are chosen to 
correspond approximately to quasi-circular (low eccen- 
tricity) inspiral. For equal-mass evolutions performed 
with the CCATIE code, parameters for quasi-circular or- 
bit were determined by minimizing an effective potential 
for the binary [23|, [67], [68[ . For the unequal-mass simu- 
lations performed with the BAM code [15|, initial mo- 
menta were specified by the 3PN-accurate quasi-circular 
formula given in Section VII of [ll[. For the longer 
unequal-mass simulations performed with higher-order 
spatial finite-difference methods [3lJ and used for verifi- 
cation, the initial momenta were taken from a PN pre- 
scription that takes radiation reaction into account to 
reduce the initial eccentricity to below e sa 10 -3 j38j . 

The Einstein equations are solved numerically with 
standard finite-difference techniques. Spatial derivatives 
are calculated at fourth- or sixth-order accuracy, and the 
time evolution is performed with a fourth-order Rungc- 
Kutta integration. Mesh refinement is used to achieve 
high resolution around the punctures and low resolutions 
far from the black holes, allowing the outer boundary to 
be placed very far (at least > 300M) from the sources. 
Full details of the numerical methods used in the two 
codes are given in [U for BAM and [H for CCATIE. 

In the wave-zone, sufficiently far away from the source, 
the spacetime metric can be accurately described as a 
perturbation of a flat background metric. Let h a b de- 
note the metric perturbation where a, b denote space- 
time indices, and t be the time coordinate used in the 
numerical simulation to foliate the spacetime by spatial 
slices. Working in the transverse-traceless (TT) gauge, 
all the information about the radiative degrees of free- 
dom is contained in the spatial part hij of h a b, where i,j 
denote spatial indices. Let us use a coordinate system 
(x,y,z) on a spatial slice so that the z-axis is parallel 
to the total angular momentum of the binary system at 
the starting time. Let 1 be the inclination angle from the 
z-axis, and let <p be the phase angle and r the radial dis- 
tance coordinates so that (r, t, <fi) are standard spherical 
coordinates in the wave-zone. 

The radiative degrees of freedom in h a b can be written 
in terms of two polarizations h + and h x : 



= h + (e + )ij + h x (e x )ij, 



(2.2) 



where e +iX arc the basis tensors for transverse-traceless 
tensors in the wave frame 



In our numerical simulations, the gravitational waves 
are extracted by two distinct methods. The first one uses 
the Newman- Penrose Weyl tensor component \E , 4 [69|, [7(| 
which is a measure of the outgoing transverse gravita- 
tional radiation in an asymptotically flat spacetime. In 
the wave-zone it can be_written in terms of the complex 
strain h = h + — \h, 



as 



lim 



71 



t' 

dt' I dt"* 4 



(2.4) 



An alternative method for wave extraction determines 
the waveform via gauge-invariant perturbations of a 
background Schwarzschild spacetime, via the Zerilli- 
Moncrief formalism (see [72| for a review). In terms of 
the even (Q^J and odd [Q^ m ) parity master functions, 
the gravitational wave strain amplitude is then given by 



1 



1 Im 



(2.5) 

Results from the BAM code have used the Weyl tensor 
component \&4 and Eq. (|2.4p . with the implementation 
described in [ll|. While the CCATIE code computes 
waveforms with both methods, the AEI-CCT waveforms 
used here were computed using the perturbative extrac- 
tion and Eq. (|2.5p . Beyond an appropriate extraction 
radius (that is, in the wave-zone), the two methods for 
determining h are found to agree very well for moving- 
puncture black-hole evolutions of the type considered 
here 

It is useful to discuss gravitational radiation fields 
in terms of spin- weighted s = — 2 spherical harmonics 
Y £ s m , which represent symmetric tracefree 2-tensors on a 
sphere, and in this paper we will only consider the dom- 
inant £ = 2, m = ±2 modes (see [28| for the higher £ 
contribution in the unequal-mass case), with basis func- 
tions 
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-^(l + cos t ) 2 e 2i *. 
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(2.6) 



Our 'input' numerical relativity waveforms thus corre- 
spond to the projections 

h tm = (Yf m 2 , h) = d<f> r h Y^f sin 9 d9 , (2.7) 
Jo Jo 

of the complex strain h, where the bar denotes complex 
conjugation. In the cases considered here, we have equa- 



(e + )ij = iilj-(f>i(f>j , and (e x )y = ti4>j + ij<i>i . (2.3) torial symmetry so that h 22 = h 2 - 2 , and 



Here I and </> are the unit vectors in the t and <p direc- 
tions, respectively, and the wave propagates in the radial 
direction. 



Ithz^ ((! + cost ) 2 h aa (t) + (1 - cos/,) 2 h aa (*)) 
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(2.8) 
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In this paper, we assume that the binary is optimally- 
oriented, so that l = 0. Thus 



h W =4 \j — h 22 {t) « 0.6308 h 22 (t). (2.9) 



III. POST-NEWTONIAN WAVEFORMS AND 
INTRODUCTION TO DATA-ANALYSIS 
CONCEPTS 

In this Section we will introduce notation that will be 
used later in the paper and describe briefly the main 
data-analysis techniques currently used in gravitational 
wave astronomy. 



A. Restricted post-Newtonian waveforms 

We use the restricted PN waveform at mass- 
quadrupole order, which has a phase equal to twice the 
orbital phase up to highest available order in the adia- 
batic approximation, and amplitude accurate up to lead- 
ing order. The corresponding h is given by 



r 



(1 + cos L) 2 e~' lip{t) + (1 - cos ife 1 ^ 

(3.1) 

where M = M x + M 2 is the total mass, 77 = M 1 M 2 /M 2 is 
the symmetric mass ratio, r is the observation radius, t 
is the inclination angle; the quantity v(t) is an expansion 
parameter, defined by v = (Mip/2) 1 / 3 with ip{t) equal to 
twice the adiabatic orbital phase. The waveform seen by 
the detector is given by 

IvI 

s(t) = 4?]— Av 2 {t)cos[<p{t) +ip ], (3.2) 
r 

where, for short-lived signals (i.e., with duration much 
shorter than the earth rotation time, as well as de- 
phasing time scale due to Doppler shifts induced by earth 
motion and rotation) , A and tpo are numerical constants 
depending on the relative position and orientation of the 
source relative to the detector, as well as the antenna 
pattern functions of the detector. In PN theory, the adi- 
abatic phase tp(t) is determined by the following ordinary 
differential equations (also called the phasing formula): 



dip 
~dt 



2v 3 
If' 



dv 
dt 



ME'(v) 



(3.3) 



In these expressions, E'(v) = dE(v)/dv where E(y) is 
the binding energy (per unit mass) of the system, and 
T(v) is the GW luminosity. E(v) and F{v) are com- 
puted as post-Newtonian expansions in terms of v [731 ]. 
Currently, the binding energy function E(v) has been 
calculated to v 6 (3PN) accuracy by a variety of meth- 
ods 0, [zi [il, tz3, m H Ea m- Thc nux Ac- 
tion !F{v), on the other hand, has been calculated to 
v 7 (3.5PN) accuracy [HI, [HJ up to now only by the 



multipolar-post-Minkowskian method and matching to 
a post- Newtonian source [73| . 

Thc inspiralling phase is usually pushed up to the 
point where the adiabatic evolution of circular orbits 
breaks down due to the lack of further stable circular or- 
bits. In the test-mass limit, the last (or innermost) stable 
circular orbit (ISCO) can be computed exactly (at 6M in 
Schwarzschild coordinates). For comparable-mass bina- 
ries, on the other hand, the ISCO cannot always arise 
unambiguously from PN theories. In adiabatic mod- 
els, the maximum-binding-encrgy condition (referred to 
as MECO, or the maximum binding energy circular or- 
bit, (83[) can be used in place of the ISCO. This condition 
is reached when the derivative of the orbital binding en- 
ergy with respect to orbital frequency vanishes. As a 
consequence, in this paper, the waveforms are evolved 
in time up to MECO: E'(v) = 0. It may be noted that 
the ISCO and MECO may not be physically meaningful 
beyond the test-mass limit, but they make convenient 
cutoff criteria. The appropriate region of validity of PN 
waveforms can only be determined by comparison with 
fully general relativistic results, such as the numerical 
simulations that we discussed earlier. 

Given E(v) and !F{v), one can construct different, but 
equivalent in terms of accuracy, approximations to the 
phasing by choosing to retain the involved functions or 
to re-expand them. Indeed, the different PN models 
which describe the GW signal from inspiralling binaries 
agree with each other in the early stages of inspiral; but 
start to deviate in the late inspiral. Thc classification 
and explicit form of various models is nicely summarized 
in [46| . In this paper we use PN waveforms obtained 
by numerically solving Eqs. (|3.3p . called the TaylorTl 
approximant, to construct thc 'hybrid waveforms' (see 
Section UVU). 



B. Introduction to matched filtering 

Since we can model the signal reasonably well, it is 
natural to employ matched filtering (which is thc optimal 
detection strategy for a signal of known shape in the sta- 
tionary Gaussian noise) to search for the gravitational- 
wave signal. Suppose the detector's data x(t) contains 
noise n(t), and possible signal s(t), i.e., x(t) = n(t)+s(t). 
Assuming n to be stationary Gaussian noise, it is con- 
venient to work in the Fourier domain, because the sta- 
tistical property of the noise is completely characterized 
by its power spectral density S n (f), which is given by 
(here we use a single- sided spectrum) 

(h(f)h*(.n)^^s n (f)s(f~.n, (3.4) 

where n(f) is the Fourier Transform of n(t) 

poo 

n{f)= / n^e-^dt, (3.5) 



6 



and (. . .) denotes taking the expectation value. Based on 
the detector noise spectrum, we introduce a Hermitian 
inner product: 



(g\h) = 2 



~9*(f)Kf) + ~g(f)h*(f) 

Sn(f) 



df. 



(3.6) 



For the data x with known signal s, the optimal detection 
statistic is given by applying a template h with the same 
shape as s, or h = as: 



Popt = (x\h). 



(3.7) 



The detectability of the signal is then determined by the 
SNR offt.pt, 



(s\h) 



N ^((h\n)(n\h)) 



(s\s. 



1/2 



(3.8) 



h—as 



(Note that the SNR does not depend on the overall nor- 
malization of h.) In case the template h is not exactly 
of the same shape as s, the SNR will be reduced to 



(3.9) 



where M. < 1 is the match of the template to the signal, 
given by 



M[s,h] 



(s\h) 



(h\h) 



(§\h). 



(3.10) 



and where a hat denotes a normalized waveform. For 
more details, wc refer the reader to Rcf. 1841. 



C. Template banks, effectualness and faithfulness 

We now consider the more realistic problem of 
attempting to detect a family of waveforms s(9), 
parametrized by a vector of physical parameters 9 € 0, 
using a family of templates /i(A) parametrized by a 
vector of parameters A <S A. We first introduce the 
concepts of physical template bank and phenomenolog- 
ical template bank. Roughly speaking, physical tem- 
plate banks are constructed from well-motivated phys- 
ical models (e.g., approximation up to a certain or- 
der) [85| , while phenomenological banks are constructed 
in an ad-hoc manner to mimic the desired physical sig- 
nals with high accuracy. For physical banks, the vectors 
9 and A consists of the same set of physical parame- 
ters, while for phenomenological banks, the vector A usu- 
ally contains phenomenological parameters, which can be 
larger or smaller in number than the physical parame- 
ters. Two phenomenological template families [H, H3] 
are used currently in the search for BH binaries in LIGO 
data [H, [8i| . They each represent a different motiva- 
tion for introducing phenomenological banks: (i) when 
we have uncertainty in the signal model, we can produce 



a template bank with larger detection efficiency by in- 
troducing extra (phenomenological) parameters (BCV1, 
(86| ) so that dim(A) > dim(8); (ii) when the true sig- 
nal depends on too many parameters and is too diffi- 
cult to search over, it is sometimes possible to come up 
with a model with fewer (phenomenological) parameters 
|dhn(A) < dim(9)) and still high fitting factors (BCV2, 

The detection efficiency of a template bank towards 
a specific signal s(0) can be measured by the thresh- 
old SNR above which the detection probability exceeds 
a certain minimum (usually 50%), while the false-alarm 
probability is kept below a certain maximum (usually 1% 
for one-year data). The threshold value depends (loga- 
rithmically, in the case of Gaussian noise) on the num- 
ber of statistically independent templates, and (inverse- 
proportionally) on the fitting factor (FF) [9(j : 

FF[h;0] = m&xM[s(9),h(\)} = M[s(9),h(X max )} . 

rnr - l ( 3 - 11 ) 

A bank with high FF is said to be effectual [J,l9l|]. Typ- 
ically, we require that the total mismatch between the 
template and true signal (including the effects of both 
the fitting factor and the discreteness of the template 
bank) to not exceed 3%. We shall see that this require- 
ment is easily met by our template bank. 

It is natural to associate every point 9 in the physical 
space O with the best matched point A max € A. This 
leads to a mapping P : 9 i— ► A defined by 



P(0) = A n 



(3.12) 



This mapping will play a key role in the construction 
of our template bank. We will assume the mapping P 
to be single- valued, i.e., given a target signal, the best- 
matched template is unique. We depict this mapping 
schematically in the left panel of Fig. [TJ 

For a physical template bank with 8 and A the same 
set of parameters (which we use 8 to denote), it is most 
convenient to identify the best-match parameter max 
as the estimation of the original parameter 8. In general 
this will lead to a systematic bias 



A0 = max - 8 = P{8) - 8 . 



(3.13) 



A bank with a small bias (as defined above) is said to be 
faithful [l,[9l|]. 

However, if we assume no uncertainty in the true wave- 
forms (thereby excluding the case of BCV1), then as long 
as P is invertible, a non-faithful physical or phenomeno- 
logical bank can always be converted into a faithful bank 
by the re-paramctrization 



'faithful 



(9) = hoP{9) 



(3.14) 



where we have used the standard notation h o P(9) := 
h(P(9)). In other words, each template A in the image 
set of physical signals P(O) is labeled by physical pa- 
rameters 9 = P _1 (A). For this reason, we require P to 
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FIG. 1: Construction of the phenomenological template bank: (i) mapping physical signals (solid curve) into a sub-manifold 
(dashed curve, with example templates marked by dots) of a larger-dimensional template bank (curved surface) , (ii) obtaining 
a lower-dimensional phenomenological bank with the same number of parameters as physical parameters, through interpolation 
(solid curve on the curved surface, with example templates marked by triangles), and (iii) Estimating the bias of the lower- 
dimensional interpolated bank by mapping physical signals into the bank (with images of example signals marked by dots). 



be invertible. It is quite conceivable that for physical 
banks, P should be invertible, if the physical bank does 
not fail to describe the true waveforms too dramatically 
(and of course assuming the true waveform does contain 
independent information about the physical parameters 
6). In this way, all reasonable physical banks can be made 
faithful. 

By contrast, if for some phenomenological bank (e.g., 
BCV2 if we only take into account the intrinsic pa- 
rameters of the bank), P is a many-to-one map, with 
P(0i) = P(0 2 ) for some 6 X ^ 2 . Then for a physi- 
cal signal with parameter Q\ , the template bank ft.f a ithful 
would achieve the same best match at both 9\ and 2 , 
making physical parameter determination non-unique. 
In this case, we can simply keep using the phenomeno- 
logical bank /i(A); once a detection is made with A max , 
the a set of parameters P _1 (A max ) would be the best 
knowledge we have about the physical parameters of the 
source. (In practice, statistical uncertainty also applies 
to A max .) 

IV. A PHENOMENOLOGICAL TEMPLATE 
FAMILY FOR BLACK-HOLE COALESCENCE 
WAVEFORMS 

A. Strategy for constructing the phenomenological 

bank 

In our situation, since it is expensive to generate the 
entire physical bank of templates using numerical simula- 
tions, we first construct a highly effectual 10-dimcnsional 
phenomenological bank (motivated by the format of PN 
waveforms), with effcctualness confirmed by computing 
its FF with a relatively small number of 'target signals'. 
Since we are considering only non-spinning black holes, 
the physical parameter space is the set of all masses 



and symmetric mass ratios (M, rf) that we wish to con- 
sider. As we shall see shortly, for our case the phe- 
nomenological parameter space A is a 10-dimensional 
space. Our templates will be denoted by 

h{\) = /hod(A) . (4.1) 

According to the discussion above [Eqs. (|3.12p - (|3.14p ]. if 
the mapping P : Q h- > A can indeed be obtained and in- 
verted, then a faithful two-dimensional (2D) phenomeno- 
logical bank can be constructed as 

/4d Mul (0) = ^ioD o P(0) . (4.2) 

However, if our aim was to know P exactly, then in 
principle we would have to calculate accurate numerical 
waveforms for every (M, rj) and to calculate the corre- 
sponding A in each case. This is obviously not practical, 
and we shall instead compute P at a few chosen points 
in and interpolate to obtain an approximation to P. 
The detailed steps are as follows: 

i. While confirming effectualness of the ten-dimensional 

(10D) bank, we simultaneously obtain TV (a num- 
ber manageable in terms of computational costs) 
data points for the mapping P, 

^L = P(0 {n) ), n=l,2,...,N (4.3) 

which gives discrete points on the 2D manifold 
P(6). This is depicted by the left panel of Fig. ED 

ii. Using these discrete points, we perform a smooth in- 

terpolation of P denoted by P nt . The form of P n t 
is motivated by PN waveforms, but with expansion 
coefficients determined by interpolation: 

Pnt(0) = A int . (4.4) 




FIG. 2: NR waveforms (thick/red), the 'best-matched' 3.5PN waveforms (dashed/black), and the hybrid waveforms 
(thin/green) from three binary systems. The top panel corresponds to rj = 0.25 NR waveform produced by the AEI-CCT 
group. The second, third and fourth panels, respectively, correspond to 77 = 0.25, 0.22 and 0.19 NR waveforms from produced 
by the Jena group. In each case, the matching region is —750 < t/M < —550 and we plot the real part of the complex strain 
(the '+' polarization). 



This gives us a 2D phenomenological bank, 

h 2 D{e) = h WD o P int (9) . (4.5) 

This is depicted by the middle panel of Fig. [T] Due 
to the discrete choice of target waveforms, the con- 
strained form of Pint, and numerical errors (in the 
target waveforms as well as in searching for best- 
fit parameters), the interpolation will have errors, 
even at the sample points. This means the 2D bank 
will have slightly lower effectualness than the 10D 
bank. 



that there will be a new mapping P^d which maps 
the physical parameters to the best fit param- 
eters of this 2D bank. We therefore find the 
best-matched parameters A^ x ,, therefore obtain- 
ing discrete samples of the mapping P 2 d : 

\tl^2D(0 (n) ), (4.6) 
yielding a systematic bias of 

A0(")=P i - 1 (A^ x ,)-0W. (4.7) 



iii. We re-test the effectualness of this 2D bank. Note 



This is depicted in the right panel of Fig. [TJ 
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In this paper, we construct the 2D template bank 
h>2T>(0) and estimate the systematic bias AO^ in the 
estimation of parameters 6, as described above. But. it 
is also possible to construct an interpolation P2D int from 
the data points of P2D so that we can construct a fully 
faithful (no systematic bias) bank (up to interpolation 
error) 



faithful 
l 2D 



(8) = h WD O P int O P 2 Dint(0)- 



(4.8) 



B. Constructing the 'target signals' 



The ultimate aim of this work is to create a family 
of analytical waveforms that are very close to the gravi- 
tational waveforms produced by coalescing binary black 
holes. As a first step, we need to construct a set of 'tar- 
get signals' containing all the three (inspiral, merger and 
ring down) stages of the binary black hole coalescence. 
Although numerical relativity, in principle, is able to pro- 
duce gravitational waveforms containing all these stages, 
the numerical simulations are heavily constrained by 
their high computational cost. It is therefore necessary, 
at the present time, to use results from post-Newtonian 
theory to extend the waveforms obtained from numerical 
relativity. 

We produce a set of 'hybrid waveforms' by matching 
the PN and NR waveforms in an overlapping time in- 
terval t\ < t < t-x- The obvious assumption involved in 
this procedure is that such an overlapping region exists 
and that in it both approaches yield the correct wave- 
forms. These hybrid waveforms are assumed to be the 
target signals that we want to detect in the data of GW 
detectors. 

The NR and PN waveforms are given by Eq. (|2.8|) 
and Eq. (|3.ip , respectively (with 1 = 0). The (com- 
plex) time-domain waveform h(t, /x) from a particular 
system is parametrized by a set of 'extrinsic parame- 
ters' fj, = {tfo,to}, where tpo is the initial phase and 
to is the start time of the waveform. We match the PN 
waveforms h PN (£, n) and the NR waveforms h NR (£, fi) by 
minimizing the integrated squared absolute difference, 5, 
between the two waveforms, i.e., 




sr In In it 

h (t, fi) — a h (i, fi) 



dt. 



(4.9) 



hybrid waveforms are then produced by combining the 
'best-matched' PN waveforms and the NR waveforms in 
the following way: 

h h y b (t, M ) = a r(t) h NR (t, M ) + (1 - r(t)) h PN (t, /i ) 

(4.10) 

where fiQ and ao denote the values of fi and a for which 
S is minimized, and r is a weighting function, defined as 



(4.11) 



In this paper we use two families of hybrid wave- 
forms. Both are produced by matching 3.5 PN Tay- 
lorTl waveforms with NR waveforms. The first set is 
constructed by using long (> 10 inspiral cycles) NR 
waveforms. This include equal-mass (77 = 0.25) NR 
waveforms produced by the AEI-CCT group using their 
CCATIE code employing fourth-order finite differencing 
to compute spatial derivatives, and equal and unequal- 
mass (?] = 0.19,0.22,0.25, or Mi/M 2 = 1,2,3) wave- 
forms produced by the Jena group using their BAM 
code employing sixth-order finite differencing and PN- 
motivated initial-data parameters. The second set of 
hybrid waveforms is constructed by using NR waveforms 
produced by the Jena group using their BAM code em- 
ploying fourth-order finite differencing. These are short 
waveforms (~ 4 inspiral cycles) densely covering a wide 
parameter range (0.16 < 77 < 0.25). We use the second 
set of hybrid waveforms to construct the phenomenolog- 
ical family and to test its efficiency in detecting signals 
from black hole coalescences, and use the first set of hy- 
brid waveforms (which are closer to the actual signals) 
to verify our results. 

The former family of hybrid waveforms is shown in 
Fig. O The NR waveforms from three different sim- 
ulations (77 = 0.25,0.22,0.19) done by AEI and Jena 
groups are matched with 3.5PN inspiral waveforms over 
the matching region -750 < t/M < -550. The hybrid 
waveforms are constructed by combining the above as 
per Eq. (|4~TUj) and Eq. ||CTT]) . 



The minimization is carried out over the extrinsic pa- 
rameters fi of the PN waveform and an amplitude scaling 
factor a, while keeping the 'intrinsic parameters' (M and 
7/) of both the PN and NR waveforms the same 2 . The 



2 Here the amplitude scaling factor a is introduced because of 
two reasons, (i) The short NR waveforms used to construct the 
phenomenological template family (see the following discussion 
in this Section) were extracted at a finite extraction radius. This 
introduces some error in the amplitude of the NR waveforms, (ii) 



Since the 'long and accurate' NR waveforms (see the following 
discussion) arc extrapolated to an infinite extraction radius, we 
expect the amplitude of these waveforms to be correct within 
numerical accuracy of the simulations. But, it turns out that the 
restricted PN waveform has an amplitude which is inconsistent 
with the NR waveform by roughly constant factor 6 ± 2% in the 
frequency range we consider here |3fj . For simplicity, we take the 
amplitude of the restricted PN waveform as the amplitude scale 
for the hybrid waveforms. It should be noted that, since we use 
normalised templates, the errors that we introduce by this (< 
10%) do not affect the fitting factors or the detection statistic. 
But the horizon distance that we estimate in Section IIV Fl can 
have an error up to 10% due to this choice. 
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The robustness of the matching procedure can be 
tested by computing the overlaps between hybrid wave- 
forms constructed with different matching regions. If the 
overlaps are very high, this can be taken as an indication 
of the robustness of the matching procedure. A p relimi- 
nary illustration of this can be found in Ref. I92j . and a 
more detailed discussion will be presented in [931 ] - 

Fig. [3] shows the hybrid waveforms of different mass- 
ratios in the Fourier domain. In particular, the panel 
on the left shows the amplitude of the waveforms in the 
Fourier domain, while the panel on the right shows the 
phase. These waveforms are constructed by matching 
3.5PN waveforms with the long NR waveforms produced 
by the Jena group. In the next section, we will try to 
parametrize these Fourier domain waveforms in terms of 
a set of phcnomcnological parameters. 



C. Parametrizing the hybrid waveforms 

We propose a phenomenological parametrization to 
the hybrid waveforms in the Fourier domain. Template 
waveforms in the Fourier domain are of particular pref- 
erence because, (i) a search employing Fourier domain 
templates is computationally inexpensive compared to 
one using time domain templates (ii) parametrization of 
the hybrid waveforms is easier in the Fourier domain. 

We take our motivation from the restricted post- 
Newtonian approximation to model the amplitude of the 
inspiral stage of the hybrid waveform, i.e., the amplitude 
is approximated to leading order as a power law /" 7 ^ 6 
in terms of the Fourier frequency / (as follows straight 
from adding leading order radiation reaction to Newto- 
nian dynamics). The amplitude of the merger stage is 
empirically approximated as a power law /~ 2 / 3 (consis- 
tent with the observation of [27]). while the amplitude 
of the ring down stage is known to be a Lorentzian func- 
tion around the quasi-normal mode ring down frequency. 
Similarly, we take our motivation from the stationary 
phase approximation (see, for example, [§H) of the in- 
spiral waveform to write the Fourier domain phase of the 
hybrid waveform as a series expansion in powers of /. As 
we shall see later, this provides an excellent approxima- 
tion of the phase of the hybrid waveform. 



1. Phenomenological waveforms 

We write our phenomenological waveform in the 
Fourier domain as 



as 



(4.12) 



where A e ff(f) is the amplitude of the waveform in the 
frequency domain, which we choose to write in terms of 
a set of 'amplitude parameters ' a = {/ me rg, /ring, o - , /cut} 



Aes(f) = C 



(/ I /merg) 
(// /merg) 



-7/6 



-2/3 



if / < fn 



if /"merg — f ^ fri 



> w £(/, /ring, o) if /ring < / < /cut 

(4.13) 

where / cu t is the cutoff frequency of the template and 
/merg is the frequency at which the power-law changes 
from /~ 7 / 6 to /~ 2 / 3 (as noted previously in [27} for the 
equal-mass case) . C is a numerical constant whose value 
depends on the relative orientations of the interferometer 
and the binary orbit as well as the physical parameters 
of the binary (see below). Also, in the above expression, 



£(/, /ring, Cr) 



1 



a 



2* J (/ - /ring) 2 + a 2 /4 ' 



(4.14) 



represents a Lorentzian function of width a centered 
around / r i n g- The normalization constant w is chosen in 
such a way that ^4 c ft(/) is continuous across the 'transi- 
tion' frequency /ring, i.e., 



TTC- / /ring 
2 \ /merg 



-2/3 



(4.15) 



Taking our motivation from the stationary-phase ap- 
proximation of the gravitational-wave phase, we write 
the effective phase ^ e s(f) as an expansion in powers of 
/, 



*eff(/) = 27r/to + W)+X)^fc/ 



(fc-5)/3 



(4.16) 



where to is the time of arrival, tpo is the frequency- 
domain phase offset, and = {?/>o, "02, ^3, "04, "06, "07 } 
are the 'phase parameters', that is the set of phenomeno- 
logical parameters describing the phase of the waveform. 

The numerical constant C in Eq. (|4. 13|) can be deter- 
mined by comparing the amplitude of the phenomeno- 
logical waveforms with that of the restricted post- 
Newtonian waveforms in the frequency domain. 

In the restricted post-Newtonian approximation, the 
Fourier transform of the gravitational signal from an 
optimally-oriented binary located at an effective distance 
d can be written as in Eq. 1|B1[) . We expect that in the 
inspiral stage (/ < / me r g ) of our phenomenological wave- 
forms the amplitude will be equal to that of the post- 
Newtonian waveforms as given in Eq. (|B1[) . Thus, in 
the case of an optimally-oriented binary, the numerical 
constant C can be computed as 



C = 



Af^/meVg 6 f5 V 



drr 2 / 3 



24 



1/2 



(4.17) 



This 'physical' scaling will be useful when we estimate 
the sensitivity of a search using this template family (see 
Section IIVFI and Appendix [5]) . 
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FIG. 3: Fourier domain magnitude (left) and phase (right) of the (normalized) hybrid waveforms. The constant phase term 
and the term linear in time (and frequency) have already been subtracted from the phase. Symmetric mass-ratio 77 of each 
waveform is shown in the legends. These waveforms are constructed by matching 3.5PN waveforms with the long NR waveforms 
produced by the Jena group. 



0.9985 
0.998 




0.16 0.18 0.2 0.22 0.24 

FIG. 4: Fitting factors of the hybrid waveforms with the 
phenomenological waveform family. Horizontal axis shows 
the symmetric mass ratio of the binary. Fitting factors are 
calculated assuming a white noise spectrum, and hence are 
independent of the mass of the binary. 



We now compute the fitting factors of the hybrid 
waveforms with the family of phenomenological wave- 
forms by maximizing the overlaps over all the param- 
eters, i.e., {at,ip,tpo,to} of the phenomenological wave- 
forms. While doing this, we also find the parameters, 
oc max and «/) max , of the 'best-matched' phenomenologi- 



cal waveforms. This calculation is described in detail in 
Appendix IA1 

We first take a few (seven) hybrid waveforms coarsely 
spaced in the parameter range 0.16 < rj < 0.25, and 
compute the fitting factors and the best-matched phe- 
nomenological parameters, assuming a white-noise spec- 
trum for the detector noise. We use these samples in the 
parameter space to construct the interpolated template 
bank (see next subsection). We then test the effectu- 
alness and faithfulness of the template bank using all 
(~ 30) hybrid waveforms finely spaced in the parameter 
space. 

The fitting factors are shown in Fig. 2J It is quite 
apparent that the fitting factors arc always greater 
than 0.99, thus underlining the effectiveness of the 
phenomenological waveforms in reproducing the hybrid 
ones. Also, as an example, in Fig. [5l we plot the hybrid 
waveforms from 77 = 0.25 binary in Fourier domain along 
with the 'best-matched' phenomenological waveform. 



2. From phenomenological to physical parameters 

It is possible to parametrize the phenomenological 
waveforms having the largest overlaps with the hybrid 
waveforms in terms of the physical parameters of the 
hybrid waveforms. In Fig. [SI we plot the amplitude 
parameters a max of the best-matched phenomenologi- 
cal waveforms against the physical parameters of the 
binary. Similarly, the phase parameters i/3 max of the 
best-matched phenomenological waveforms are plotted 
against the physical parameters of the binary in Fig. [7J 
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FIG. 5: Hybrid waveforms (solid lines) in the frequency domain, and the 'best-matched' phenomenological waveforms (dashed 
lines). The left panel shows the Fourier domain magnitude, while the right one shows the phase. In the hybrid waveforms, 
the constant phase term and the term linear in time (and frequency) have already been subtracted from the phase. In the 
phenomenological waveforms, to and ipo (see Ea. (|4,16|) ) have been chosen to minimize the phase difference between the hybrid 
and phenomenological waveforms. These waveforms correspond to a binary with r\ = 0.25, and are constructed from the 
'short' NR waveforms produced by the Jena group (see Section TlVB \ . The 'dip' in the left panel at Mf ~ 2 x 10~ 2 is due to 
the small eccentricity present in the first few cycles of the NR waveform. All waveforms are normalized assuming a flat noise 
spectral density. 



It can be seen that a max and i/> max can be written as 
quadratic polynomials in terms of the physical parame- 
ters (M and rf) of the hybrid waveforms as: 



*j int 



1>, 



k int 



Q j V 2 + M + Cj 
nM 

XkV 2 +ykV + z k 
7 ? (7rM)( 5 - fc )/ 3 



(4.18) 



where dj, b 3 , 



Cj, j = 0...3 and x k , Vk, Zk, k 
0. 2, 3, 4, 6, 7 are the coefficients of the quadratic poly- 



nomials used to fit the data given in Figs. [6] and [7l These 
coefficients are listed in Tables III and ITT1 It may be noted 
at this point that Figs. [5] and [7] correspond to the map- 
ping P : — ► Amax that we have introduced in Sec- 
tion IIV Ai and Eq. (|4.18[) to the interpolation P mt of 
P. 



Using the empirical relations given in Eq. (|4.18p . we 
can rewrite the effective amplitude and phase of the 
waveforms in terms of M and rj as: 



a ri 2 +b ri+c 



A eS (f) = C I 



iMf 

o»?+co 



-7/6 



-2/3 



a-aV +bo?;+co 
ttM 



if /< 

■r a ri 2 +bai]+co f - axif+biri+ci 
11 TrAf — J 



*eff(/) 



ya a ri 2 +b (l 

w r ( f aiy 2 +bir)+ci a.2ri 2 +t>2r)+C2 \ aii-j 2 +bn]+ci ^ r 

1 7 

27r/t + <po + - Y,( x k V 2 + Vk V + Zk) (ttM/) 



tM 



< 



a3V~ +b3V+C3 

7TM 



\(fc-5)/3 



(4.19) 



fc=0 



where the constant C is given by Eq. (|4.17p . We use this family of parametrized waveforms to create a two- 
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FIG. 6: Best-matched amplitude parameters a max in terms of 
the physical parameters of the binary (assuming white noise 
spectrum). The horizontal axis shows the symmetric mass- 
ratio of the binary. Quadratic polynomial fits oti nt to the 
data points are also shown. 



dimensional template bank of non-spinning waveforms. 
This template family can be seen as a two-dimensional 
sub-manifold (parametrized by M and 77) embedded in 
a higher dimensional manifold (of the phenomenological 
waveforms). 

The polynomial coefficients in the Table HT1 are indeed 
significantly different from those predicted by stationary 
phase approximation of the PN inspiral phase in the fre- 
quency domain. There are two reasons for that: The 

I 



where f 



where / 



where /o = 215 Hz. 

Faithfulness is computed by maximizing the overlaps 
over the extrinsic parameters to and <po only, which can 
be done trivially |96| . Effectualness is computed by max- 
imizing both intrinsic and extrinsic parameters of the 



first one is that our re-parametrization is optimized for 
the mass range where all three phases (inspiral, merger 
and ring down) are contributing significantly. The sec- 
ond reason is the residual eccentricity present in the nu- 
merical waveforms. Change in the relative significance of 
different PN terms reflects attempt to match the slightly 
eccentric waveform with circular. When more accurate 
(less eccentric) numerical waveforms become available in 
future, the re-parametrization given in Eq. (|4.18p can be 
optimized for a wider mass range. An example of this 
can be seen in Ref. [92| . 

D. Effectualness and faithfulness 

In order to measure the accuracy of our parametrized 
templates we compute their overlap with the 'target sig- 
nals' (the hybrid waveform). To check the faithfulness 
of our phenomenological templates, we compute their 
overlap with the target signal maximizing it over the ex- 
trinsic parameters (time-of-arrival and the initial phase) . 
We assess the effectualness of the parametrized wave- 
forms by computing fitting factors with the target sig- 
nals (computing the overlap maximized over both extrin- 
sic and intrinsic parameters). Faithfulness is a measure 
of how good the template waveform is in both detect- 
ing a signal and estimating its parameters. However, 
effectualness is aimed at finding whether or not an ap- 
proximate template model is good enough in detecting 
a signal without reference to its use in estimating the 
parameters. 

We compute the effectualness and the faithfulness of 
the template family for three different noise spectra. The 
one-sided noise power spectral density (PSD) of the Ini- 
tial LIGO detector is given in terms of a dimcnsionlcss 
frequency x = f/f by [H] 



(4.20) 



(4.21) 



(4.22) 



binary. The maximization over the intrinsic parameters 
is performed with the aid of the Nelder-Mead downhill 
simplex algorithm [97| . 

The effectualness of the template waveforms with the 
hybrid waveforms is plotted in Fig. [5] for three different 



Sh(f) = 9 x 1CT 46 [(4A9x)- 56 + 0.16x~ 4 - 52 + 0.52 + 0.32a; 2 ] 
= 150 Hz; while the same for Virgo reads [95| 

S h (f) = 10.2 x 10~ 46 [(7.87x)- 4 ' 8 + 6/Ux- 1 + 1 + x 2 ] , 
= 500 Hz. For Advanced LIGO H, 

'l-x 2 +x 4 /2^ 



S h (f) = 10 



-4!) 



x- 414 - 5x- 2 + III 



/l — X + X / Z\ 

I l + a; 2 /2 ) 
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FIG. 7: Best-matched phase parameters i/> m ax in terms of the physical parameters of the binary (assuming white noise 
spectrum). The horizontal axis shows the symmetric mass-ratio of the binary. Quadratic polynomial fits ipi nt to the data 
points are also shown. 



noise spectral densities. The corresponding faithfulness 
is plotted in Fig. [9j It is evident that, having both values 
always greater than 0.99, the proposed template family 
is both effectual and faithful. 

We also calculate the systematic bias in the estima- 
tion of parameters while maximizing the overlaps over 
the intrinsic parameters of the binary. The bias in the 
estimation of the parameters is defined in Eq. ()4.7j) . 

The percentage biases in estimating the total mass 
AI, mass ratio 77, and chirp mass M c = Mrf/ b of the 
binary are plotted in Figs. [TUl HU and \Y2\ respectively. 
This preliminary investigation suggests that the bias in 
the estimation of M and rj using the proposed template 
family is < 3%, while the same in estimating M c is < 6%. 



Parameter 


afe 


b k 


Cfe 


/merg 


2.9740x10" 


1 


4.4810x10" 


2 


9.5560x10" 


2 


/ring 


5.9411x10" 


1 


8.9794x10" 


-2 


1.9111x10" 


1 


a 


5.0801x10" 


r 


7.7515x10" 


-2 


2.2369x10" 


■2 


/cut 


8.4845x10" 


1 


1.2848x10" 


-r 


2.7299x10" 


1 



TABLE I: Polynomial coefficients of the best-matched ampli- 
tude parameters. The first column lists the amplitude param- 
eters CKint- Eq. (|4.18|l shows how these parameters are related 
to the coefficients ak,bk,Ck. 



Parameter 




Vk 




V>o 


1.7516xl0 _1 


7.9483 xlO -2 


-7.2390xl0 -2 


1p2 


-5.1571X10 1 


-1.7595X10 1 


1.3253X10 1 


lp3 


6.5866xl0 2 


1.7803 xlO 2 


-1.5972xl0 2 


V>4 


-3.9031 xlO 3 


-7.7493xl0 2 


8.8195xl0 2 


^6 


-2.4874x10* 


-1.4892xl0 3 


4.4588 xlO 3 


1p7 


2.5196xl0 4 


3.3970xl0 2 


-3.9573xl0 3 



TABLE II: Polynomial coefficients of the best-matched phase 
parameters. The first column lists the phase parameters i/'int- 
Eq. (|4.18[) shows how these parameters are related to the co- 
efficients x k ,y k ,z k . 



E. Verification of the results using more accurate 
hybrid waveforms 

As we have discussed in Section llV B[ the hybrid wave- 
forms used for constructing the template waveforms are 
produced by matching rather short (~ 4 inspiral cycles) 
NR waveforms with PN waveforms. We have also pro- 
duced a few hybrid waveforms by matching PN wave- 
forms with long (> 10 inspiral cycles) and highly accu- 
rate (sixth-order finite differencing and low eccentricity) 
NR waveforms. This set of hybrid waveforms (which are 
closer to the 'actual signals') can be used to verify the 
efficacy of the template waveforms in reproducing these 
more accurate signals. 
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FIG. 8: Fitting factor of the (two-dimensional) template family computed using three different noise spectra. The panel in 
the left correspond to the Initial LIGO noise PSD, the one in the middle to the Virgo noise PSD and the one in the right 
to the Advanced LIGO noise PSD. The horizontal axis represents the symmetric mass ratio 77 of the binary and the legends 
display the total mass M (in units of Mq). 
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FIG. 9: Same as in Fig. [8] except that the plots show the faithfulness of the (two-dimensional) template family. 



Fig.[T5lshows the fitting factors of the two-dimensional 
template family with the 'more accurate' hybrid wave- 
forms. The fitting factors arc computed, as before, using 
the Initial LIGO (left), Virgo (middle) and Advanced 
LIGO (right) noise spectra. The high fitting factors (al- 
though smaller than the same obtained in the previous 
Section) with the hybrid waveforms once again underline 
the efficacy of the template waveforms in reproducing 
the hybrid ones. It is indeed expected that the tem- 
plate family will have better overlaps with the hybrid 
waveforms described in the previous Section (those con- 
structed from 'short' NR waveforms), because the poly- 
nomial coefficients given in Tables HI and ILTl are optimized 



for these hybrid waveforms. When more 'long and ac- 
curate' NR waveforms are available in the future, the 
polynomial coefficients given in the Tables can be op- 
timized for the corresponding family of 'more accurate' 
hybrid waveforms. In any case, since the fitting factors 
are already very high, we don't expect any significant 
improvements. 
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M M M 

FIG. 10: Bias in the estimation of M. Horizontal axis reports the total mass M (in units of M@) and vertical axis reports the 
symmetric mass ratio r\ of the binary. Colors in the plot corresponds to the percentage bias, \AM\/M x 100. The left panel 
corresponds to the Initial LIGO noise PSD, the middle panel to the Virgo noise PSD and the right panel to the Advanced 
LIGO noise PSD. 
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FIG. 11: Same as in Fig. 1101 except that the plots show the percentage bias \ A-q\/r] x 100 in the estimation of rj. 



F. The astrophysical range and comparison with 
other searches 



The template family proposed in this paper can be 
used for coherently searching for all the three stages 
(inspiral, merger, and ring-down) of the binary black 
hole coalescence, thus making this potentially more sen- 
sitive than searches which look at the three stages sep- 
arately. Fig. [M] compares the sensitivity of the searches 
using different template families. What is plotted here 
are the distances at which an optimally-oriented, equal- 
mass binary would produce an optimal SNR of 8 at the 
Initial LIGO (left plot), Virgo (middle plot) and Ad- 
vanced LIGO (right plot) noise spectra. In each plot, 
the thin solid (blue) line corresponds to a search using 



PN templates truncated at the innermost stable circu- 
lar orbit (ISCO) of the Schwarzschild geometry having 
the same mass as the total mass M of the binary; the 
dashed (purple) line to a search using ring-down tem- 
plates [98( ; the dot-dashed (black) line to a search us- 
ing effective one body Q waveform templates truncated 
at the light ring of the corresponding Schwarzschild ge- 
ometry, and the solid line to a search using all three 
stages of the binary coalescence using the template bank 
proposed here. The computation is described in detail 
in Appendix [Bj The horizontal axis reports the total- 
mass of the binary, while the vertical axis the distance 
in Mpc. It is quite evident that, for a substantial range 
of total mass (100 < M/M Q < 300 for Initial LIGO, 
200 < M/Mq < 500 for Virgo, 150 < M/M Q < 400 for 
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FIG. 12: Same as in Fig. 1101 except that the plots show the percentage bias \ AM C \/M C x 100 in the estimation of M c = Mr] 3 / 5 
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FIG. 13: Fitting factor of the (two-dimensional) template family with 'more accurate' hybrid waveforms (see Section [IV El) . 
The overlaps are computed using three different noise spectra. The panel in the left correspond to the Initial LIGO PSD, the 
one in the middle to the Virgo noise PSD and the one in the right to the Advanced LIGO noise PSD. The horizontal axis 
represents the total mass M (in units of M©) and the legends display the symmetric mass ratio r\ of the binary. 



Advanced LIGO), the 'coherent search' using the new 
template family is significantly more sensitive than any 
other search considered here. 



However, while this looks promising, we repeat here 
the caveats emphasized in [30| : It is important to treat 
Fig. [14] as only a preliminary assessment; fitting factors 
are not the only consideration for a practical search strat- 
egy. It is also very important to consider issues which 
arise when dealing with real data. For example, false 
alarms produced by noise artifacts might well determine 
the true sensitivity of the search, and these artifacts will 
inevitably be present in real data. This is however be- 
yond the scope of the present work, and further investi- 
gation is required before we can properly assess the effi- 
cacy of our phenomenological template bank in real-life 
searches. 



V. SUMMARY AND OUTLOOK 



Making use of the recent results from numerical rel- 
ativity we have proposed a phenomenological waveform 
family which can model the inspiral, merger and ring- 
down stages of the coalescence of non-spinning binary 
black holes in quasi-circular orbits. We first constructed 
a set of hybrid waveforms by matching the NR wave- 
forms with analytical PN waveforms. Then, we con- 
structed analytical phenomenological waveforms which 
approximated the hybrid waveforms. The family of phe- 
nomenological waveforms that we propose was found to 
have fitting factors larger than 0.99 with the hybrid 
waveforms. We have also shown how this phenomeno- 
logical waveform family can be parametrized solely in 
terms of the physical parameters (M and rj) of the bi- 
nary, so that the template bank is, in the end, two di- 
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FIG. 14: Distance to optimally located and oriented equal-mass binaries which can produce an optimal SNR of 8 at the Initial 
LIGO (left plot), Virgo (middle plot) and Advanced LIGO (right plot) noise spectra. Horizontal axis reports the total mass of 
the binary (in units of M©) and vertical axis reports the distance in Mpc. In each plot, the thin solid (blue) line corresponds 
to a search using standard PN templates truncated at ISCO, the dotted-dashed (black) line to a search using effective one 
body waveform templates truncated at the light ring, the the dashed (purple) line to a search using ring-down templates, and 
the thick solid (red) line to a search using the template family proposed in this paper. The ring down horizon distance is 
computed assuming that e = 0.7% of the black hole mass is radiated in the ring down stage, while the Kerr parameter a — 0.69 
is known from the numerical simulation. Since the value of e has some amount of uncertainty in it, we have also included the 
shaded region in the plot corresponding to 0.18% < e < 2.7%. 



mensional 3 . This two dimensional template family can 
be explicitly expressed in terms of the physical param- 
eters of the binary. We have estimated the 'closeness' 
of this two-dimensional template family with the fam- 
ily of hybrid waveforms in the detection band of three 
ground-based GW detectors, namely Initial LIGO, Virgo 
and Advanced LIGO. We have estimated the effectual- 
ness (larger overlaps with the target signals for the pur- 
pose of detection) and faithfulness (smaller biases in the 
estimation of the parameters of the target signals) of the 
template family. Having both types of overlap always 
greater than 0.99, the two dimensional template family 
is found to be both effectual and faithful in the detection 
band of these ground-based detectors. 

This phcnomenological waveform family can be used 
to densely cover the parameter space, avoiding the com- 
putational cost of generating numerical waveforms at 
every grid point in the parameter space. We have com- 
pared the sensitivity of a search using this template fam- 
ily with other searches. For a substantial mass-range, the 
search using all three stages of the binary black hole co- 
alescence was found to be significantly more sensitive 
than any other template-based searches considered in 
this paper. This might enable us to do a more sensitive 
search for intermediate-mass black holes using ground- 



3 It may be noted that, the mapping from the phenomcnologi- 
cal to physical parameters might not be unique in the case of 
spinning binaries, because of the degeneracies of different spin 
configurations. 



based GW detectors. 

A number of practical issues need to be addressed be- 
fore we can employ this template family in an actual 
search for GW signatures. The first issue will be how to 
construct a bank of templates sufficiently densely spaced 
in the parameter space so that the loss in the event rate 
because of the mismatch between the signal and tem- 
plate is restricted to an acceptable amount (say, 10%). 
The explicit frequency domain parametrization of the 
proposed template family makes it easier to adopt the 
formalism proposed by Owen [99| in laying down the 
templates using a metric in the parameter space. Work 
is ongoing to compare the metric formalism adopted to 
the proposed template family and other ways of laying 
out t he te mplates, for example a 'stochastic' template 
bank 100]. Also, this explicit parametrization makes it 
easier to employ add ition al signal-based vetoes, such as 
the 'chi-square test' |l0l| . This will also be explored in 
a forthcoming work. 

Since this template bank is also a faithful representa- 
tion of the target signals considered, we expect that, for 
a certain mass-range, a search which coherently includes 
all three stages of the binary coalescence will bring about 
remarkable improvement in the estimation of parameters 
of the binary. This may be especially important for LISA 
data analysis in estimating the parameters of supermas- 
sive black hole binaries. This is also being explored in 
an ongoing work [47j . 

It is worth pointing out that the family of target sig- 
nals (the hybrid waveforms) that we have considered 
in this paper is not unique. One can construct alter- 
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nate families of hybrid waveforms by matching PN wave- 
forms computed using different approximations with NR 
waveforms. Also, owing to the differences in initial data 
and accuracy of numerical techniques, the NR waveforms 
from different simulations can also be slightly different. 
Thus, the coefficients listed in Tables Q] and [XT] have a 
unique meaning only related to this particular family 
of target waveforms. But we expect that the general 
parametrization that we propose in this paper will hold 
for the whole family of non-spinning black hole coales- 
cence waveforms from quasi-circular inspiral. As we have 
mentioned in the Introduction, the purpose of this paper 
is to explicitly prescribe a general procedure to construct 
interpolated template banks using parametrized wave- 
forms which mimics actual signals from binary black hole 
coalescence (as predicted by numerical relativity and an- 
alytical methods). 

Nevertheless, it may be noted that most of the PN 
waveforms constructed using different approximations 
are known to be very close to each other (see, for ex- 
ample, [ill). Also, we expect that NR waveforms from 
different simulations will converge as the accuracy of n u- 
merical simulations improves (see, for example, jl02j |). 
Thus, since different families of PN and NR waveforms, 
which are the 'ingredients' for constructing our target 
signals, are very close to each other, we expect that the 
phcnomenological waveform family proposed in this pa- 
per, in its present form, will be sufficiently close to other 
families of target signals for the purpose of detecting 
these signals. As a preliminary illustration of this, we 
have computed the fitting factors of the template wave- 
forms with a different family of hybrid waveforms (con- 
structed from longer and more accurate NR waveforms), 
and have shown that the overlaps are indeed very high. 
This will be explored in detail in a forthcoming work. 

Also, we remind the reader that this paper consider 
only the leading harmonic of the GW signal {1 = 1, m = 
±2). We expect that the contribution from the higher 
harmonics become important for high mass ratios, which 
will be investigated in a forthcoming work. 
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APPENDIX A: CALCULATION OF THE 
FITTING FACTORS 

In order to find the fitting factor of our phenomeno- 
logical bank to a hybrid waveform, as well as the 
best-matched parameters (a max ,i/) max ), we need to per- 
form a maximization of the overlap A4(a, ip) in a 12- 
dimensional space, which seems a challenging task at 
first sight, especially due to the oscillatory nature of the 
dependence of A4(ot, ip) on the components of How- 
ever, due to the very high fitting factor, as well as the 
linear dependence of ty c s (*/>', /) onip, we have been able 
to design an analytic approximation to A4(a, tp) that is 
highly accurate and can be maximized over ip analyti- 
cally. In describing this approximation, we also include 
ipo and to in tp, forming an 8-dimcnsional vector. 

For a target hybrid waveform 



and a phcnomenological template 

u(f) = A eB ( a ;f)e i *««&->n 



(Al) 



(A2) 



the overlap M(a,ip) can be broken into a product of 
two terms, 



M(a,ip) = Ma(oc) Mp{a,^) 



(A3) 



with 



W5 irw, ( A4) 



and 



M P (a, *) = 1 r ^(";/W)cos[A*(/)] d/ 



where 



Sh(f) 



A*(/) = *(/)-*eff(V;/). 



(A5) 



(A6) 



In the above expressions, the normalization constants a 
and b are defined by 



a 2 = 



A 2 {f) 
Sh(f) 



df 



A 2 cS (*-J) 
Sh(f) 



df, 



(A7) 
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and 



1. Search using post-Newtonian templates 



b = 



A cS {a-f)A{f) 
S h (f) 



(A8) 



If the phase difference A\E'(/) is small, we can approx- 
imate cos »1- A^ 2 /2, and rewrite A4p as 



Mp « M'p = I 



26 



A eff (a;/)A(/)[A*(/)]< 



S h (f) 



df- 

(A9) 

Since ^S c s(ip; f) is a linear function in xp, minimizing A4' P 
becomes a least-square fit with a weighting function 



Kf) 



A cS (a;f)A(f) 
Shif) 



(A10) 



More specifically, writing 1 J r e ff(V'i /) as m Eq. (|4.16|l . i.e., 

* eff (V;/) = E^/ (5 ^ )/3 ' ( An ) 



we have 



1-Mp = - [ipAi> T -2Bip T + D] , (A12) 

where we have defined a matrix A, a vector B and a 
scalar constant D, such that 

1 Z" 00 

^ = iy / (io -^ )/3 M/)d/, 



'0 

= / (5 " 3)/3 ^(/)/i(/)d/. 
D = -J *2 (/)M/)d/ . 



The maximum of A4' P is then equal to 



M'p 



reached at 



1 — - to — B A -1 Bl 



BA- 1 . 



(A13) 



(A14) 



(A15) 



As a consequence, for each a, we are able to maximize 
Mp(a,tp), and hence A4(a,i/>), over xp analytically. 
The original 12-dimcnsional maximization is then con- 
verted to a 4-dimensional maximization, only over the 
amplitude parameters, on which the overlap depends in 
a non-oscillatory way. 



In the restricted post-Newtonian approximation, the 
Fourier transform of the gravitational signal from an 
optimally-oriented binary located at an effective distance 
d can be written in the following way: 



Kf) 



dn 2 / 3 



5 // 
24 



1/2 



-7/6 e i[27i-/t - 



•¥>o+W)-V4] 



(Bl) 

where M is the total mass, r\ is the symmetric mass 
ratio, to is the time of arrival and ipo is the initial phase. 
The phase ip(f) is computed using the stationary phase 
approximation. 

The optimal SNR in detecting a known signal h buried 
in the noise is given by 



Popt — 2 



df 



Kff 



s h {f) 



1/2 



(B2) 



where Sh(f) is the one-sided PSD of the noise. The 
optimal SNR in detecting the signal given in Eq. (|Bl[) 
can thus be computed as: 



Popt 



drr 2 / 3 



5 i] 
IT 



1/2 



df 



f 



-7/3 



Shif) 



1/2 



(B3) 



where /i ow is the low-frequency cutoff of the detector 
noise and / U pp is upper frequency cutoff of the template 
waveform. The effective distance to a binary which can 
produce an optimal SNR p opt can be computed by in- 
verting the above equation. 

The standard post-Newtonian waveforms are trun- 
cated at / upp = /isco, where / IS co = (6 3/2 7rA/) _1 is 
the GW frequency corresponding to the innermost sta- 
ble circular orbit (ISCO) of the Schwarzschild geome- 
try with mass equal to the total mass M of the bi- 
nary. The effective one body (EOB) waveforms are trun- 
cated at / upp = / L r, where / L r = (3 3/2 7rA/) _1 is the 
GW frequency corresponding to the light ring of the 
Schwarzschild geometry with mass M. Both of these 
quantities are computed assuming the test particle limit. 
It may be noted that, for the EOB waveforms, an analyt- 
ical Fourier domain representation is not available. They 
cannot be expressed in the form given in Eq. (|Bl|) . But 
for the purpose of the estimation of the horizon distance, 
these formulas give a reasonable approximation. 



2. Search using ring down templates 



APPENDIX B: COMPUTING THE HORIZON 
DISTANCE 

Here we describe how we compute the horizon dis- 
tance of different searches discussed in Section HVFl An 
alternative way of computing the horizon distance can 
be found in Ref. [§1. 



The ring down portion of the GW signal from a coa- 
lescing binary, considering only the fundamental quasi- 
normal mode, cor respo nds to a damped sinusoid. This 
can be written as [103| 



ftring(t) = A ring exp 



7t/qnr(£ - tp) 

Q 



x cos (-27r/ QNR (i - t ) + <fo) , (B4) 



21 



where A r - mg is the amplitude, to is the start time of the 
ring down, 930 the initial phase, M is the mass of final 
black hole, / QNn and Q are the central frequency and the 
quality factor of the ringing. For the fundamental mode, 
a good fit to the frequency / QNR and quality factor Q, 
within an accuracy of 5%, is given by 



QNR 



[1- 0.63(1 -a) 3/10 ] 



1 



2irM ' 



Q » 2(1 -a)" 9 / 20 , 



(B5) 
(B6) 



where aM 2 is th e spi n angular momentum, and a is the 
Kerr parameter [1031 ] . 

To compute the optimal SNR in dete cting this signal 
present in the data, we proceed as in [104| . assuming 
that for t < to, faring (i) is identical to t > to except for 
the sign in the exponential, and dividing by a correcting 
factor of s/2 in amplitude to compensate for the doubling 
of power: 



faring (0 



exp 



V2 

x cos (-2nf QNR (t - to) + ip ) ■ (B7) 



Its Fourier transform then becomes 

L . ( f\ _ ^rmg J QNR i27T/t [ T 

rmsUj " V2nQ U 2 + 4(/-/ QNR ) 2 



.9 2 + 4(/ + / QNR )V ' 



(B8) 



where g = f QNR /Q. 

In general, it is not easy to estimate A r i ng , or the two 
polarization amplitudes; they depend upon the detailed 
evolution of the merger epoch, as well as variables such 
as the orientatio n of the final m erged remnant. A reason- 
able hypothesis pMllM[l07| is that their ratio follows 
the ratio of the inspiral polarization amplitudes. With 
this hypothesis, the overall amplitude of the signal from 
an optimally located and oriented binary, requiring that 



the ring down radiate some fraction e of the system's 
total mass, becomes 



^opt = /5eM 2 
Iing V ^ d y/ Mf Qm QF(Q) 



(B9) 



where F{Q) = 1 
source. The optimal SNR p can now be computed as 



24 q 2 and d is the distance to the 



Popt 



/up 



df 



\hn 



"I 1/2 



ingl 



fh 



s h (f) 



(BIO) 



where /i ow and / upp are the lower and upper cutoff fre- 
quencies of the detector noise. As in the previous case, 
the horizon distance can be computed by inverting this 
equation. 

3. Search using the template family proposed in 
this paper 

The phcnomcnological waveforms in the frequency do- 
main are given in Eqs. (|4.12H4.17|) . The optimal SNR in 
detecting this signal can be computed as: 



Popt 



M 5 / 6 / mc 7 / g 6 fh 

d7T 2 /3 



/ring 



df 



(///. 



mcrg; 



1/2 



-4/3 



/mcrg I J m( , rg ) 7 / 3 



df 



Sh(f) 



Sh(f) 



/out /: 2 (/,/rin g ,a) 
/ rinE S h (f) 



-,1/2 



(Bll) 



where £(/, / r i, lg , ct) is defined in Eq. Q4.14[) . and 
/mcrg, /ring, /cut and u are given by Eq. (|4.18p . 

This equation can be inverted to calculate the effec- 
tive distance to the optimally-oriented binary which can 
produce an optimal SNR p Q pt- 
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